rm(list=ls())
set.seed(123456)

png(filename = "05testplot.png", width=6, height=6.5, units="in", res=400)

crit <- qt(0.975, df=100)

x <- seq(from=-1, to=8, by=0.01)
null.plot <- dt(x, df=100)
plot(null.plot ~ x, type="n", ylab="probability density", xlab="estimated t statistic",
     ylim=c(0, 0.55))

x.shade <- seq(from=crit, to=8, by=0.01)
y.shade <- dt(x.shade, df=100, ncp=3)
polygon(c(min(x.shade),x.shade), c(min(y.shade), y.shade), col=gray(0.8, alpha=1), border=NA)

y.shade <- dt(x.shade, df=100)
polygon(c(min(x.shade),x.shade), c(min(y.shade), y.shade), col=gray(0.5, alpha=1), border=NA)

null.plot <- dt(x, df=100)
lines(null.plot ~ x)

alt.plot <- dt(x, df=100, ncp=3)
lines(alt.plot ~ x, lwd=2)

abline(v = crit, lty=2)

arrows(5.5, 0.3, 3, 0.15, length=0.1)
power.value <- round( pt(crit, df=100, ncp=3, lower.tail = F) + pt(-crit, df=100, ncp=3, lower.tail = T), 4)
text(6, 0.325, labels=(paste( c("power = ",  power.value), collapse = " "  ) ) )

arrows(5.25, 0.15, 2.5, 0.025, length=0.1)
size.value <- round( pt(crit, df=100, lower.tail = F) + pt(-crit, df=100, lower.tail = T), 4)
text(6.5, 0.175, labels=(paste( c("size = ", size.value), collapse = " "  ) ) )

legend("topright", lty=c(1,1,2), lwd=c(1,2,1), pt.cex=1, cex=0.8,
       legend=c("t density under null", "t density, ncp = 3", "critical-t, alpha = 0.05"))

dev.off()



rm(list=ls())
set.seed(123456)

png(filename = "005testplot.png", width=6, height=6.5, units="in", res=400)

crit <- qt(0.9975, df=100)

x <- seq(from=-1, to=8, by=0.01)
null.plot <- dt(x, df=100)
plot(null.plot ~ x, type="n", ylab="probability density", xlab="estimated t-statistic",
     ylim=c(0, 0.55))

x.shade <- seq(from=crit, to=8, by=0.01)
y.shade <- dt(x.shade, df=100, ncp=3)
polygon(c(min(x.shade),x.shade), c(min(y.shade), y.shade), col=gray(0.8, alpha=1), border=NA)

y.shade <- dt(x.shade, df=100)
polygon(c(min(x.shade),x.shade), c(min(y.shade), y.shade), col=gray(0.5, alpha=1), border=NA)

null.plot <- dt(x, df=100)
lines(null.plot ~ x)

alt.plot <- dt(x, df=100, ncp=3)
lines(alt.plot ~ x, lwd=2)

abline(v = crit, lty=2)

arrows(5.5, 0.3, 3.75, 0.15, length=0.1)
power.value <- round( pt(crit, df=100, ncp=3, lower.tail = F) + pt(-crit, df=100, ncp=3, lower.tail = T), 4)
text(6, 0.325, labels=(paste( c("power = ",  power.value), collapse = " "  ) ) )

arrows(5.25, 0.15, 3, 0.025, length=0.1)
size.value <- round( pt(crit, df=100, lower.tail = F) + pt(-crit, df=100, lower.tail = T), 4)
text(6.5, 0.175, labels=(paste( c("size = ", size.value), collapse = " "  ) ) )

legend("topright", lty=c(1,1,2), lwd=c(1,2,1), pt.cex=1, cex=0.8,
       legend=c("t density under null", "t density, ncp = 3", "critical-t, alpha = 0.005"))

dev.off()







rm(list=ls())
set.seed(123456)

png(filename = "power-comparison.png", width=6, height=6.5, units="in", res=400)

size <- seq(from=0.5, to=0.9975, by=.0001)
crit <- qt(size, df=100)
power <- pt(crit, df=100, ncp=3, lower.tail = F) + pt(-1*crit, df=100, ncp=3, lower.tail = T)
plot(power ~ size, type="l", ylim=c(0,1))

size <- seq(from=0.5, to=0.9975, by=.0001)
crit <- qt(size, df=100)
power <- pt(crit, df=100, ncp=2, lower.tail = F) + pt(-1*crit, df=100, ncp=2, lower.tail = T)
lines(power ~ size, type="l", lty=2)

size <- seq(from=0.5, to=0.9975, by=.0001)
crit <- qt(size, df=100)
power <- pt(crit, df=100, ncp=1, lower.tail = F) + pt(-1*crit, df=100, ncp=1, lower.tail = T)
lines(power ~ size, type="l", lwd=2, col="gray")

legend("bottomleft", lty=c(1,2,1), lwd=c(1,1,2), col=c("black", "black", "gray"),
       legend=c("beta = 3", "beta = 2", "beta = 1"))

dev.off()



png(filename = "multiple-test-power.png", width=6, height=6.5, units="in", res=400)

coef <- seq(from=0.01, to=6, by=0.01)
crit <- qt(0.9975, df=100)
power <- pt(crit, df=100, ncp=coef, lower.tail = F) + pt(-1*crit, df=100, ncp=coef, lower.tail = T)
crit.2 <- qt( 1-(0.005^(1/2))/2, df=100)
power.2 <- pt(crit.2, df=100, ncp=coef, lower.tail = F) + pt(-1*crit.2, df=100, ncp=coef, lower.tail = T)
power.2 <- power.2^2
crit.3 <- qt( 1-(0.005^(1/3))/2, df=100)
power.3 <- pt(crit.3, df=100, ncp=coef, lower.tail = F) + pt(-1*crit.3, df=100, ncp=coef, lower.tail = T)
power.3 <- power.3^3

plot(power.3 ~ coef, type="l", lwd=2, col="gray", ylim=c(0,1),
     ylab = "power", xlab = expression(paste("true relationship ", beta)))
lines(power ~ coef)
lines(power.2 ~ coef, lty=2)

legend("bottomright", lty=c(1,2,1), lwd=c(1,1,2),
       col=c("black", "black", "gray"), 
       legend = c("one test", "two tests", "three tests"))

dev.off()









# make a plot of bias associated with lowered statistical significance
# for a specific relationship, beta = 3, sd_hat-beta = 1, df = 100
# some code used from: https://www.r-bloggers.com/setting-graph-margins-in-r-using-the-par-function-and-lots-of-cow-milk/


rm(list=ls())
set.seed(123456)

png(filename = "pub-bias-single.png", width=6, height=6.5, units="in", res=400)

alpha.v <- seq(from=0.005, to=0.1, by=0.001)
bias.v <- rep(NA, length(alpha.v))
crit.v <- qt( 1-(alpha.v/2), df=100 )

samples <- rt(1e7, ncp=3, df=100)

pb <- txtProgressBar(min=1, max=length(crit.v), initial=0, style=3)
for(i in 1:length(crit.v)){
  
  setTxtProgressBar(pb, i)
  bias.v[i] <- mean( (samples[ which(samples > crit.v[i] ) ]  / 3) - 1 )
  
}
close(pb)

bias.v <- 100*bias.v

plot(bias.v ~ alpha.v, type="l", xlim=c(0, 0.1), ylim=c(0, 30),
     ylab = "estimated publication bias (as % of coefficient)",
     xlab = "alpha for statistical significance test")

dev.off()

# generate a plot for an overall prior distribution including a spike
# at beta = 0

rm(list=ls())
set.seed(123456)

png(filename = "pub-bias-dist.png", width=6, height=6.5, units="in", res=400)

gen.norm.samp <- function(se.b, df.b, pr.false, n.draws, bound, seed){
  
  if(seed!="off"){
    set.seed(seed)
  }
  b.true <- rnorm(n.draws, mean=0, sd=bound)
  b.false <- rep(0, n.draws)
  choice <- ifelse(runif(n.draws)<(1-pr.false), 1, 0)
  b.test <- ifelse(choice==1, b.true, b.false)
  normalizer <- mean(abs(b.test))
  
  b.test.est <- b.test + se.b*rt(n.draws, df=df.b)
  bias <- b.test.est - b.test
  
  mult <- sign(b.test)
  mult <- ifelse(mult==0, sign(bias), mult)
  bias.out <- (mult*bias)/normalizer
  
  b.t <- abs(b.test.est/se.b)
  
  list.out <- list(b.true = b.test, b.est = b.test.est, b.bias = bias.out, b.t = b.t)
  return(list.out)
  
}


alpha.v <- seq(from=0.005, to=0.1, by=0.001)
bias.v <- rep(NA, length(alpha.v))
crit.v <- qt( 1-(alpha.v/2), df=100 )

samples.list <- gen.norm.samp(se.b=1, df.b=100, pr.false=0.5, n.draws=1e7, bound=3, seed=123456)

pb <- txtProgressBar(min=1, max=length(crit.v), initial=0, style=3)
for(i in 1:length(crit.v)){
  
  setTxtProgressBar(pb, i)
  sig <- which(samples.list$b.t > crit.v[i])
  bias.v[i] <- mean(samples.list$b.bias[sig])
  
}
close(pb)

par(mar=c(5,6,4,2)+0.1,mgp=c(4,1,0))

bias.v <- 100*bias.v

# add line to the plot
plot(bias.v ~ alpha.v, type="l", xlim=c(0, 0.1), ylim=c(40, 50),
     ylab = c("estimated publication bias","(as % of mean relationship in prior)"),
     xlab = " ")
mtext(text = "alpha for statistical significance test", side=1, line=3)

dev.off()

